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ABSTRACT 

We investigate numerically the propagation of density waves excited by a 
low-mass planet in a protoplanetary disk in the nonlinear regime, using 2D lo- 
cal shearing box simulations with the grid-based code Athena at high spatial 
resolution (256 grid points per scale height h). The nonlinear evolution results 
in the wave steepening into a shock, causing damping and angular momentum 
transfer to the disk. On long timescales this leads to spatial redistribution of the 
disk density, causing migration feedback and potentially resulting in gap opening. 
Previous numerical studies concentrated on exploring these secondary phenom- 
ena as probes of the nonlinear wave evolution. Here we focus on exploring the 
evolution of the basic wave properties, such as its density profile evolution, shock 
formation, post-shock wave behavior, and provide comparison with analytical 
theory. The generation of potential vorticity at the shock is computed analyti- 
cally and is subsequently verified by simulations and used to pinpoint the shock 
location. We confirm the theoretical relation between the shocking length and 
the planet mass (including the effect of the equation of state), and the post-shock 
decay of the angular momentum flux carried by the wave. The post-shock evolu- 
tion of the wave profile is explored, and we quantitatively confirm its convergence 
to the theoretically expected N-wave shape. The accuracy of various numerical 
algorithms used to compute the nonlinear wave evolution is also investigated: we 
find that higher order spatial reconstruction and high resolution are crucial for 
capturing the shock formation correctly. 

Subject headings: Planet-disk interactions, protostellar disks. Hydrodynamics, 
Methods: numerical, Planets and satellites: formation 



^ Sloan Fellow 



-2- 



INTRODUCTION 



In recent years the study of disk-planet interaction in protoplanetar y disks, which aims at 
expla i ning; the enormous dive rsity of the observed exoplanetary systems (jPapaloizou &: Terquem 
20061 : iPapaloizou et al.l 120071 ). has received much attention. Specifi cally, the discoveries o f 
Jovian pla nets and brown dwa r fs residing at both ext remely tight (IMayor &: Queloa 119951 ) 
and wide ( IChauvin et al.l l2004l : iNeuhauser et al.l |2005| ) orbits challenge the current models 
of planet formation, and inspire studies of planetary migration driven by disk-planet inter- 
action. 

A planet embedded in a protoplanetary disk launches spiral density waves in the disk, 
which propagate away from the planet. Gravitational interaction between these nonaxisym- 
metric density waves and the planet exchanges angular momentum between them, driv- 
ing type I planetary migration for a low mass planet, which can not open a gap in the 
disk. The timescale for t ype I migration is u s ually very short ( < 10^ vr for nlanfits mnrp 



massive than a few Mm, Tanaka et al. 2002 



D'Aneelo & Lubow 


20081) 


xlO^ yr. 


Hartmann et al. 



19981 ). which makes planetary survival a major puzzle. Sever al solutions to this problem 
have been proposed, involving magnetic fields (iTergueml 120031). turbulent density fluctua- 
tions c aused by magnetorotation al instability (MRI) (iNelson fc Papaloizoul 120041. density 



traps (iMenou fc Goodman 



2004r). and co-orbital torques (iPaardekooper fc Mellemal 12006 



20081 : iKlev k CridJ2008[ ). to name a 



Baruteau fc Mass"etl 20081 : Paardekooper &: Papaloizou 
few. 

However, the majority of the theoretical calculations of type I migration do not take 
i nto account the ch ange of the disk density profile caused by the planetary density waves 
( iTanaka et al.ll2002l ). In reality, density waves carry angular momentum as they travel away 
from the planet. Eventually, this angular momentum is deposited elsewhere in the disk, 
leading to the redistribution of the disk mass. For a migrating low-mass planet this density 
redistribution slightly enhances the disk surface density in front o f the planet and reduces 
the density behind it, slowing down or even halting the migration ( iHourigan fc Wardlll984j : 
Ward fc Houriganlll989t IWardlll997l : lRafikovll2002bl ). This effect is called migration feedback. 
At higher planetary masses density redistribution by planetary torques is so severe that a 
gap may form near the planetary orbit. This shifts the planetary migration from type I to 
type II, and effe ctively slows the raigratio n down, which may save planets from falling onto 
the central star ( iLin &: Papaloizoulll986al Jbl: IWardlll997| ). 



It is important to emphasize that the change of the state of the di sk by planetary density 
wave s can only be accomplished by virtue of some damping processes (iGoldreich &: Nicholson 
19891 ). which can be either linear or nonlinear. Several possible linear wave damping mecha- 
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nisms hav e been proposed, such as radiative damping (ICassen &: Woolumlll996l ) and viscous 
damping ( iTakeuchi et aLl Il996l ). However, all of these mechanisms have drawbacks and 
may not be ve ry efficient at dissipating an gular momentum carried by the planet-generated 
density waves (iGoodman &: Rafikovll200ll . hereafter GROl). 



GROl proposed another wave damping process, which results from the nonlinear wave 
evolution into a shock. This process could serve as a universal nonlinear wave damping 
mechanism working at regi mes where a l terna tives fail (specifically, it does not require a 
background disk viscosity) ( 1L arsoni 1 1 98 9l . 1 1 99 Ol . GROl). GROl analytically investigated the 
nonlinear propagation of the density wave in the local shearing-sheet approximation, and 
predicted the dependence of the shocking length Ish, the radial separation between the orbit 
of the planet and the point where t h e shock first appears, on the equation of state and the 
planet mass. Subsequently, iRafikovl (l2002al ) extended the local treatment of GROl into the 
global case to include the effects of surface density and temperature variations in the disk 
as well as the disk cylindrical geometry and nonuniform shear. 



Previously, iLin fc Papaloizoul ( 119931 ) have suggested the following condition for the gap 
opening: 



no 



(1) 



where Cg is the sound speed in the disk and flp is the angular velocity of the planet. The 
thermal mass Mth is the mass of a planet at which the Hill radius and the Bond i radius 
Rb = GM-p/ of the planet are comparable to the scale height of the gaseous disk h flRafikov 
20061 ). For a MMSN model faavashilEgSll . Minimum Mass Solar Nebulae) and = Mq 

3 / \ 3/4 



th 



12 



1 km s 



1 AU 



(2) 



where Vp is the semi-major axis of the planetary orbit. 



GROl found that planets with Mp > Mth generate density waves which are nonlinear 
from the very start, and shock as soon as they are produced. Their angular momentum is 
deposited in the immediate vicinity of the planet, within (1 ~ 2)h, leadi ng to the surface 
density evolution and gap opening close to the planet. On the other hand, IRafikovl (l2002bl ) 
demonstrated that density waves produced by lower mass planets are still able to dissipate 
efficiently further out through the nonlinear damping, even if they are only weakly nonlinear 
to begin with. In fact one should expect that given enough time, even a very low mass 
planet at a fixe d semi-major ax is should be able to open a gap in an inviscid disk. Based on 
this argument. IRafikovl (l2002bl ) explored gap opening mediated by the nonlinear dissipation 
of density waves and indeed found that the low mass planets which do not satisfy the 
condition ([1]) can be capable of stalling their migration via the migration feedback and 
opening a gap in a low viscosity disk. 
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However, most previous numerical studies have failed to capture gap opening by small 
planets (Mp < Mth). This is most likely because these simulations have studied flows with 
significant viscosity (artificial or numer ical). In a viscous disk there is another condition 
Lin fc Papaloizoulll993l : lRafikovll2002b[ ) 



th 



a a 



0.043 h 



1/2 



(3) 



that has to be fulfilled simultaneously with ([T]) for the gap to open, where the effective 
Shakura-Sunyaev a = u/hcs {v is the kinematic viscosity). It is quite likely that the levels of 
viscosity in these simulations were too high (in particular, the numerical viscosity due to the 
low spatial resolution) for the torques induced by the low mass planets to overcome viscous 
filling of the gap, i.e. the viscous criterion ^ was not fulfilled. 



Only recently has the GROl theory been confirmed in numerical work. iPaardekooper 



(120061 ) first found evidence for the distance over which density waves damp to increase with 
decreasing Mp, exactly as predicted by GROl. Subsequently, using two-dim ensional hydro- 
dynamic simulations of migration of low-mass planets in nearly inviscid disks, iLi et al.l (120091 ) 
found the migration rate to drop due to the migration feedback. They also showed that high 
disk viscosity (a > 10"'^) could wash out this effect, making fast type I migration persist. 
In addition, they confirmed the existence of a critical planet mass which eventua lly halts 
th e mig r ation, and f o und r easonable agreement with the theoretical prediction of iRafikov 



( l2002bl ). iMuto et al.l (|2010[ ) also investigated nonlinear wave evolution by 2D hydrodynam- 
ical simulations and confirmed that low mass planets (a few tenths of Mth) indeed are able 
to open gaps in an inviscid disk. They successfully detected shock formation and density 
redistribution in the disk, resulting from the nonlinear wave evolution. However, these simu- 
lations investigated only secondary consequences of the nonlinear wave evolution, such as the 
slowdown of the migration due to migration feedback and the limit on gap opening planet 
mass. 

It is obvious from this discussion that elucidating the mechanisms of nonlinear wave 
dissipation in protoplanetary disks is not only interesting in and of itself, but is also crucial 
for understanding the migration feedback and the gap opening issue, b oth of which could 



significantly increase the planetary migration timescale. Here, following iDong et al.l ( 12011 



here-after Paper I), who explored the linear excitation and propagation of planet-generated 
density waves, we continue to investigate the evolution of density waves via numerical sim- 
ulations. We now focus on the nonlinear wave evolution, and provide detailed quantitative 
comparisons with analytical results of GROl. 

The structure of this paper is as follows. In § [21 we briefiy summarize the main results of 
the nonlinear theory, and analytically study the potential vorticity generation at the shock. 
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which we use to pinpoint the shock location. The code and numerical setup of the simulations 
are introduced in § |3l We then present the main numerical results and compare them with 
theory in § HJ The effects of variation of numerical algorithms are investigated in § and 
the effect of linear damping due to viscosity is studied in § O We discuss the implications 
for realistic protoplanetary disks in § [71 and summarize main results in § |8l 



2. Nonlinear Theory of the Density Wave Propagation 



Here we briefly summarize the results of GROl for the nonlinear density wave evolution. 
For a planet of sufficiently low mass [i.e. Mp <^ Mth), the excitation and initial propaga- 
tion of the wave are linear processes, which are not affected by the nonlinearity. Far from 
the planet, the wave excitation is no longer important, while the nonlinear effects start to 



accumulate. In the shearing sheet geometry (as usual, x 



r — Tr 



and y = rp{9 — 9p) denote 



pseudo- Cartesian radial and azimuthal coordinates in a corotating system centered on the 
planet, and {u,v) represents the perturbed velocity in {x,y) plane) GROl have shown that 
under the assumption of weak nonlinearity the fully nonlinear system of fluid equation can 
be reduced to a single first-order nonlinear equation: 



dtX - xdqX = 0, 



(4) 



which is the inviscid Burger's equation ( iWhithaml Il974| ) . The dimensionless variables ap- 
pearing here are related to radius, azimuth, and density contrast as follows: 



t 
V 

X 



37/2 



211/45 
3 
2 



X 

h 

3x2 



5/2 



y 



sign(x) 



23/4(7 + 1; 
33/2 



1/2 



S - En M 



th 



(5) 
(6) 

(7) 



where h is the scale height ar id Cs is the sound sp eed of the disk, and 7 is the adiabatic 
index in the equation of state. IPaardekooperl ( 20061 ) has demonstrated that in the case of an 
isothermal disk in definition ([7]) (S — So)/So should be replaced with ln(S/So). However, 
since we are interested only in the weakly nonlinear wave evolution the assumption of (S — 
So)/So ^ 1 is always made and Eq. (JTj) holds even in the isothermal case. 

In the absence of linear damping mechanisms such as viscosity, for almost any choice of 
smooth initial conditions in Burger's Eq. (jl]), the solution will eventually become double- 
valued, which means that a shock must appear. This phenomenon is analogous to the 
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nonlinear evolution of sound waves, except that for planetary density waves (in the linear 
regime) the conservation of the wave angular momentum flux (AMF) and differential rotation 
lead to the wave amplitude increasing, and the radial wavelength decreasing with distance 
from the excitation point in the linear regime. Both effects conspire to accelerate wave 
steepening, and the formation of shocks close to the planet. GROl provided an analytical 
formula for the nonlinear shocking length l^h'- 



sh 



/7 + 1 Mp 

IT275 mT 



th 



-2/5 



(8) 



which indicates that if Mp > Mth, the shock appears very close to planet (/sh ~ h), i.e. in 
the region where the wave excitation is taking place. GROl also predicted that deep into 
the post shock regime {x ^ Zgh) the density wave profile should attain the N-wave shape 
(ILandau fc Lifshitzlll959l ). The amplitude of the N-wave scaled by a;^/^ 



^0 



1/2 



■'-'th 



and its azimuthal width w should scale with distance from the planet as: 



A oc 



W (X t 



1/2 



(9) 



(10) 



where the time-like coordinate t is defined by Eq. ([5]). Note that A oc t° (stays roughly 
constant) in the linear regime of wave evolution due to the angular momentum conservation. 

After the shock formation, the wave gradually damps out, transferring its angular mo- 
mentum to the mean fiow, and forcing the disk to evolve. The angular momentum fiux 
Fh{x) carried by the wave decays in the post-shock region as (GROl): 



Fh{x) oc \x\~^/'^ {\x\ > /sh) 



;ii) 



Where the shearing-sheet ge ometry i s assumed. In the global setting the post-shock wave 
evolution was investigated in iRafikovl (l2002al ). 



2.1. Generation of potential vorticity 

A very useful diagnostic of the density wave evolution is provided by the so-called 
potential vorticity, sometimes also called vortensity. This quantity is defined in the shearing 
box coordinates as 

^^e,.(Vxv) + 2n^ (12) 
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where v is the fluid velocity in the frame of the shearing sheet, and is the unit vector 
perpendicular to the disk plane. The background value of this quantity for v = is Co = 
fi/(2So). In an inviscid, barotropic fluid, potential vorticity is conserved everywhere, except 
across shocks. At the shock ( experiences a jump the magnitude of which depends on 
the strength of the shock and the orientation of the shock front with respect to the incoming 
flow. Thus, the potential vorticity perturbation = C ~ Co can be interpreted as the 
evidence for the appearance of shock waves and provides a useful diagnostic of the flow. 

Here we theoretically calculate the behavior of AC as a function of distance from the 
planet and planetary mass Mp based on the weakly-nonlinear theor y of GROl. We start with 
the following expression for th e jump in C at the isothermal shock ( iKe vlahani 1 1 9 9 71 : iLi et al. 
2005l : lLin fc Papaloizoulboioh : 



AC 



Cs (M2 - 1)^ dM 



/AS 



dS \ S 



_d_ 

dS 



AS\ 



(13) 



where M is the Mach number of the flow perpendicular to the shock, AS is the surface 
density jump at the shock front and 5* is the distance along the shock. To arrive at the last 



expression we used the relation 



AS/S valid in the isothermal case. 



— 1 /2 

We can write d/dS = [1 + [dysh/dxY] d/dx, where ysh{x) ~ —{3/A)h{x/hy is the 



theoretical shock position in the x — y plane. For x > h we find then 



d_ 

dS 



2h d 
3x dx 



(14) 



We now use Eq. ([7]) to relate AS/S to the jump across the shock Axsh of the function x, 
satisfying the Eq. (jl]): 



AS 



33/2 



1/2 M, 



^-Axsh{t{x)), 



(15) 



th 



S 23/4(7 + 1) 
where t{x) is given by Eq. ([5]). 

Combining Eqs. ([5]), f lT3|) - f|T5|) . and setting M ^ 1 (shock is weak) in Eq. f|T3|) . we find 

37/2 



ACix) 



213/4(7 + 1) 



dlnAxsh , ^ 



[Axshitix))Y 





h 


1/2 


/Mp\ 




X 


1 





(16) 



Note AC 7^ only for x > Igh, where the shock exists. Eq. ( [T6|) shows that potential vorticity 
generation is a steep function of Mp-. at a fixed distance x away from the planetary orbit 
AC oc M3. Thus, low mass planets are very inefficient at generating potential vorticity and 
detecting its production in numerical experiments with small planets is a non-trivial task. 



- 8 - 



Using Eq. ([8]) and setting 7 = 1 we can rewrite Eq. f|T6l) as: 
AC(x/U) ^ 1.3^ 

This particular form is natural since ^Xsh is a function of t ~ \x/lsh['^'^ (for isothermal gas) 
only, as follows from Eq. (jlj), which does not contain any free parameters. Thus, a fixed 
x/lsh corresponds to a unique value of t and /S.Xsh- 

For instance, if we choose x/lgh (and t) to correspond to the point where /S^Xsh (and A(^) 
attains its maximum value, the corresponding ^Cmax would scale with Mp as oc (Mp/Mth)"*^^^^. 
More generally, one expects {Mt\i/ M^)^^^^ to depend on x/lsh only. 

We compare these theoretical predictions with numerical results in §4.21 



Ijh 

X 



th 



iO/ o 



d In Axsfe 
(91n(x/U) 



[I\Xsh{t{x/lsh))Y 



(17) 



3. Code, Method and the Numerical setup 



To provide quantitative comparison with analytical theory in GROl and in §2.H we 
carry out 2D shearing box hydrodynamics simulations using Athe na, a grid-based code for 



astrophysical gas dyna mics using higher-order Godunov methods (iGardiner fc Stond 12005 



2008c IStone et al.ll2008l ). The basic numerical setup is the same as in paper I, and we briefly 
summarize it below. 



The implementation of the shearing box approximation in Athena is described in lStone fc Gardiner 



(120101 ). This approximation adopts a frame of reference located at radius corotating with 
the disk at orbital frequency Vtp = f2(rp), and the Cartesian coordinate system {x,y) is the 
same as in § [21 The rate of shear corresponds to a Keplerian rotation profile. The equation 
of state (EOS) we use in the simulations is the ideal gas law with 



P 



7-1 



where Eij^ is the internal energy, p is the gas pressure, and 7 = 5/3. We also use an isothermal 
EOS (effectively 7 = 1), in which case p = Sc^, where Cs is the isothermal sound speed. We 
do not include explicit viscosity except in § [6] (z.e. no linear dissipation is present in the 
system), and we do not account for the self-gravity of the disk. 

To study the accuracy and convergence of our results we have varied the numerical 
algorithms used to compute our simulations Unless noted otherwise, our standard 

simulations use the following choices (the same as in paper I): an isothermal equation of 
state, the Roe solver with third order reconstruction in characteristic variables, and the 
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corner transport upwind (CTU) unsplit integrator (IStone et al.ll2008l ). a resolution of 256 
grid points per scale height h (subsequently denoted 256//i for brevity), and fourth order 
accurate potential ^p"^^ 

$(4) = _GM ^ ^"^^^^ fl9) 

^^'^P(-p2 _^ ^2)3/2 ^^^) 

with softening length = /?./32, where p = a/o;^ + is the cylindrical radius counted from 
the position of the planet. This potential converges to Keplerian at p ^ rg as {r^/ p)^ (which 
means the fractional error is ©((r^/p)^) as r^/p — )■ 0). 

We use the following boundary conditions (BCs, the same as paper I). On x (radial) 
boundaries, we keep values of all physical variables in ghost zones fixed at their respective 
unperturbed Keplerian values {i.e. keep the ghost zones as their initial states), and the waves 
leave through the x boundaries when they reach the edge. On the y (azimuthal) boundary, we 
experimented with two BCs: the conventional outflow EC, and an inflow/outflow BC. In the 
former case the variables in the ghost zones are copied from the last actively-updated row of 
cells. In the latter case, the variables in the ghost zones are fixed at their initial values if they 
are the physical "inflow" boundaries (the regimes y < 0,x < and ?/ > 0, a; > 0), or copied 
from the last actively-updated row of cells if they are the physical "outflow" boundaries (the 
regimes y > 0,x < and y < 0,x > 0). We found that the conventional outflow y BC 
accumulates some non-zero velocity perturbation on top of the pure linear shear velocity 
profile, while our inflow/outflow y BC does not. This affects calculation of variables derived 
from the simulated velocity field, such as potential vorticity. We use the inflow/outflow y 
BC for our simulations. 

As in paper I, we check the level of numerical viscosity in our runs. We run a series of test 
simulations with otherwise identical conditions but different explicit Navier-Stokes viscosity, 
and measure the wave properties. As the explicit viscosity decreases the simulation results 
gradually converge to the one with zero explicit viscosity, which indicates the numerical 
viscosity dominates the explicit viscosity. For a typical simulation with low Mp = 2.09 x 
10~^Mth and an isothermal equation of state, the effective Shakura-Sunyaev a-parameter 
(a = vjhc^ characterizing our numerical viscosity is found to be b elow 10~^. Suc h small 



levels of viscosity are expected in dead zones of protoplanetary disks ( lGammielll996l ). where 
magnetorotational instability (MRI) may not operate effectively (?). Also note even then 
the MRI operates, MHD turbulence may not act like a Navier-Stokes viscosity. 

For the standard simulations in this work we use box size X'oh x 102/;,, thus the overall 
grid resolution in our runs is 4096 x 26112 (for the standard resolution of 256//;.). We note 
that these domains are larger than in paper I since we want to follow the nonlinear evolution 
of fluid variables deep into the post-shock regime. In a few cases with very small Mp, we 
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extend the simulation box size to 20h x 156h to fit the large 4h inside the box (see Eq. |8]) 
Our simulations are run for at least 10 and in some cases up to 100 orbital periods. 



4. Results and Discussion 

In this section, we present our numerical results on the density wave properties in the 
nonlinear regime, and make detailed quantitative comparison with the theory. 



4.1. Density Profile in the Post Shock Region 

The nonlinear evolution of the wave is depicted in Figure [IK, where we show the az- 
imuthal density profile at two pre-shock locations, the theoretically predicted shocking dis- 
tance Ish, and four post-shock positions. In accord with linear theory, before the shock the 
numerical peak amplitude of the wave stays roughly constant (we see only a slight increase 
with X expected as a result of the continuing accumulation of the torque from low azimuthal 
wavenumber harmonics by the wave). At the same time the leading edge of the density pro- 
file gradually steepens as the wave propagates away from the planet. At x ~ Ish the leading 
edge becomes vertical and the wave shocks. After that the leading edge stays vertical while 
the height of the profile A (see Fig. [1^ for an illustration of its definition) decreases due to 
the dissipation of energy and angular momentum. 

We note that the N-wave shape does not quite appear in our wave profile, because the 
second (trailing) shock does not emerge until very large distance {x ~ 7(Mth/Mp)2/5/i), which 
our simulation boxes do not cover. Nevertheless, the leading segment of the N-wave does 
appear and we use its behavior to quantitatively explore the N-wave evolution. In Figure [Th, 
we plot the dependence of A and the azimuthal width of the wave w (the y/h — 3{x/h)'^/4: 
value at the mid point of the leading edge, as indicated in Fig. [1^) on the coordinate t{x) 
defined by Eq. (|5]) in the post shock regime, as well as the theoretical scaling relations (Eq. 
[To]) with arbitrary normalization. Far in the post shock regime {t ^ 1) our numerical results 
on the N-wave behavior agree well with theory. 



Evolution of the wave profile in the nonlinear regime was first studied by iMuto et al. 



( l2010l ). who demonstrated the steepening of the wave profile and the decay of S after the 
shock formation (their Figure 6). However, they did not make quantitative comparison 
between the numerical results and theoretical scaling relations, as we do here. 
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4.2. Excitation of potential vorticity 

In Figure [2] we show a typical spatial patten of potential vorticity deviation = 
( — (fi/2So) from its background value in our simulations (only one half of the simulation 
box is shown). Because of the velocity shear, the fluid enters the box from the two inflow 
y boundaries {x > Sz y > and x < Sz y < 0) with initial background = 0, and 
maintains this value until it reaches the spiral wave. If the fluid element meets a spiral 
wave at |xcross| ^ ^sh, the shock is not crossed and the original A^ value is conserved. The 
fluid element carries it until it leaves the simulation domain through the outflow y boundary 
(x < & y > and x < k y > 0). 

On the other hand, if |xcross| ^ ^sh, the fluid element crosses the shock, potential vorticity 
conservation is broken and ( gets a kick. The amplitude of the ( jump depends on |xcross|- A 
larger IxcrossI results in a weaker ( jump because further from the planet the shock becomes 
weaker and the incidence angle of the fluid on the shock is more oblique. After the shock, 
the fluid element conserves the new value of ( until it leaves the box. Note that ( is ill 
de fined right at the shock front, which produces the feature along the spiral (also visible 



m 



Lin fc Papaloizoull2010l Fig.l), a nd the noisy structure at small \x\ is due t o the vortex 



generation in the co-orbital region (IKoller et al.l 120031 : iLin fc Papaloizoul 120101 ). Measured 



as a function of |x| at large \y\ <^ h, A( maintains the background value as small |x| until 
\x\ ~ Ish, where it first grows with x, and then gradually decreases, as the wave decays and 
the shock becomes weaker and weaker. 

In Figure |3^ we plot the maximum value of the vorticity jump A(rnax (scaled by i^/So 
to make it dimensionless) as a function of the planetary mass Mp. The value of A(fnax "was 
derived by varying x at fixed y ^ ysh(^sh) behind the shock and finding the maximum of 
A({x) at the shock in a simulation with a given Mp. Analytical calculation of A( generation 
in §2.11 based on weakly nonlinear theory of GROl predicts that A(max oc Mp^^^. However 
we find that our results are fit marginally better by the dependence A^max oc M^-^^ over 
more than 2 orders of magnitude in planetary mass. This result confirms the low efficiency 
of the potential vorticity excitation by low-mass planets. 

In Figure [3b we plot the full radial profiles of A^ produced after a single shock crossing 
for three different Mp. For better graphic representation, we smooth the ( curve by averaging 
value over /i/32 in x to reduce the noise. We display the data scaled by M^-^^ as a function 
of x/lsh- The exponent 2.95 of the mass scaling is the same as in Figure [3^ so that the 
resultant curve has a universal shape, independent of Mp. 

One can see that, as expected, A^ stays close to zero for x < Igh, i-e. prior to the 
appearance of the shock. At the shock A^ very rapidly (within radial distance < O-^hh) 
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attains its maximum value, as the shock front develops. Subsequently decreases because 
the shock rapidly weakens with increasing x: its amplitude is reduced by the dissipation 
and also the shock becomes more oblique to an incoming flow. As stated in §2] far from the 
planet shock evolves into an N-wave with x decaying as Axsh oc oc {x/lsh)~^^'^, see 

Eqs. ([5]) and (fTOj) . As a result, according to Eq. ( !T7|) . at x ^ Ish the amplitude of should 
generally decrease oc \x/l, 



sh 



-17/4 



At smaller x > Igh ^Xsh does not obey the N-wave scaling and decays with x rather 
slowly, as can be directly seen from Figure [TJj, where A acts as a proxy for /S.Xsh- Also, 
at these separations the factor in square brackets in Eq. f|T6|) varies in a non-trivial fash- 
ion. Initially dhiAxsh/d^'^t is small, see Figure [TJd, and the factor in square brackets 
(as well as the A^) is positive. However, at large separations N-wave evolution results in 
din ^Xsh/ dint — 7- —1/2 and the factor in square brackets approaches —3/2, making A^ 
negative. From our numerical calculations we flnd that A(^ cha nges sign at x k, 1.5/sh- This 
result agrees well with numerical calculations of other authors ( iLi et al.ll2005l ). 



Previously, potential vorticity generation at the shock was studied bv iKol 



( l2003l ). who focused on the co-orbital region of the protoplanet, and by iLi et al. 



er et al. 



(120051 ) 



who numerically investigated the dependence of potential vorticity on spatial resolution in 
2D inviscid disks. Both stu dies addre s sed th e flow instability caused by the potential vor- 
ticity generation. Recently, lYu et aL fcoioh explored the time evolution of the potential 
vorticity and its dependence on r^. iMuto et al.l (120 lOf ) also studied the possibility of using 
potential vorticity to identify the shock formation in an attempt to verify the theoretical 
/sh — scaling relation in GROl. However, their potential vorticity proflles were rather 
noisy preventing meaningful quantitative comparison. 



Lin fc Papaloizoul ( l2010l ) investigated the potential vorticity generation by a massive 



planet (Mp > Mth, so the wave shocks immediately after being excited). They followed 
fluid elements on horseshoe orbits, and conflrmed that the potential vorticity is generated as 
material passes through the two spiral shocks. In a global cylindrical geometry employed in 
their work, a fluid element gets a kick in the potential vorticity every time it passes a spiral 
shock, so the potential vorticity in the simulation box increases with time. This is also the 
case in our 2D local shearii ig sheet geometry when we switch the y boundary condition to 
periodic. However, while in lLin fc Papaloizoul ( l2010l ) the potential vorticity stops increasing 
and reaches a plateau after 30 — 50 orbits, in our low mass planet (Mp ^ Mth) and high 
resolution cases we do not see this saturation. In one experiment with 100 orbits, our AC, 
linearly increases with time throughout the entire simulation time. It is not clear what causes 



^In the N-wave regime one can identify /S.Xsh from Eq. (fTS]) with A defined by Eq. (|9]). 
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this difference, but we suspect that the lo w resolution (whic h intro duces larger numerical 



viscosity) might be responsible for it in the iLin &: Papaloizoul (120 lOf ) case 



4.3. Effect of Equation of State 

In linear regime, the density wave evolution does not depends on the equation of state, 
as long as the sound speed of the gas is fixed. However, the EOS has a prominent effect in the 
nonlinear regime, which results in a dependence of the shock location on 7 (Eq. [8]). We show 
in Figure H] the radial A( profiles for two simulations with otherwise identical parameters 
(including the sound speed) but with different EOS (7 = 1 and 7 = 5/3). There are two 
major differences between them. First of all, Ish for the two cases are different, with larger 
7 resulting in earlier shock. The difference between /gh in two cases is consistent with the 
theoretical prediction for these 7 (~ 10%). Second, the peak value in the 7 = 5/3 case 
is about 25% higher than in the 7 = 1 case. On the other hand, the decay of the A( profiles 
is qualitatively similar for different EOS. 



4.4. The 4h - Mp Relation 

One of the most important results of the analytical theory by GROl is the relation 
([8]) between the shocking length /gh and the mass of the planet Mp. This relation plays a 
central role in the calculation of every pr ocess dr i ven by the nonlinear evolution, such as 



the migration feedback and gap opening (IRafikovl l2002bl ) . Here we provide the numerical 



confirmation of this relation, as shown in Figure [51 

We define the shock location as the value of x at which ( reaches half of its maximum 
value at the jump. The numerical data points nicely agree with the theoretical expectation 
Ish oc Mp"^^^ for about 2.5 orders of magnitude in Mp, with deviation < 10% for most of 
the Mp range. The smallest Mp presented here (3.69 x 10^'^Mth) corresponds to ~ 4 Lunar 
mass and the largest Mp (0.667Mth) corresponds to ~ 8M® at 1 AU for a MMSN model. 
The numerical result deviates from theory at the largest Mp as expected, because the wave 
excitation and the shock formation regions are no longer separated there. At the small Mp 
end we also see a trend of deviation from theory. This is because the nonlinearity for such 
low-mass planets is so small that the linear dissipation due to numerical viscosity becomes 
non- negligible and starts to damp the wave prior to the theoretically expected location (see 
discussion in § [H]) • We expect that modeling the disk at even higher resolution than we have 
(256//i) or more accurate algorithms will resolve this problem (§[S])- 
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Previous attempts to verify Ish — Mp relation ([8]) were pioneered by iPaardekooperl ( 120061 ). 
who inferred Zgh from the width of the gap opened by the planet. They found to increase 
with decreasing Mp, in agreement with GROl theory. However, the resoluti on of his s i mula- 
tions was insufficient for quantitative verification of Eq. ([H]). Subsequently, IYu et al.l (|2010[ ) 
confirmed Ish — Mp relation using potential vorticity as means of shock detection. Their cal- 
culations spanned only about an order of magnitude in Mp and were not fully converged {e.g. 
in terms of soft ening length r^). Since their simulations were run in cylindrical geometry. 



Yu et al.l ( I2OIOI ) found Igh in the inner disk to be smaller than in the outer disk. This is to 



be expected, since global shear rate is highe r in the inner disk, which causes faster nonlinear 



evolut ion there in agreement with iRafikovl (l2002a( ). Nevertheless, the results of lYu et al. 
( 120101 ) for Ish — Mp relation are in reasonable quantitative agreement with Eq. ([8]) of this 
paper. 



4.5. The decay of Angular Momentum Flux 

The spatial patten of the decay of the angular momentum fiux carried by the density 
wave determines how the disk will respond to the presence of the protoplanets, which subse- 
quently determines the efficiency of the density feedback and the gap opening. GROl studied 
the post-shock AMP decay, and predicted the asymptotic behavior of AMP at large distance 
(Eq. ( ITT]) , also see Pigure 3 in GROl.). Numerically, we calculate the AMP as: 

Fh{x) = J:o j uvdy, (20) 

—00 

where u and v are the azimuthal and radial velocity perturbation of the fiuid with respect to 
the background shear profile. In theoretical calculation, GROl showed that Fh{x) is given 
by 

. ^JI^ ( (21) 



23/2(^+1)2^ V^th 

where the dimensionless AMP $(t) is defined as 

m= I X\v,t)dr]. (22) 



and Xy I1 cind t are defined in Eqs. ([MZ])- Here we follow their notation and measure $(t) 
from our simulations by calculating Fh{x) and using Eq. (I?I]) . 

We plot the numerically measured $(t) in Pigure [6] for three different values of Mp. Our 
numerical results agree well with the theoretical prediction at large distance after the shock. 
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which is the working range of the theory, showing that the AMF decay is indeed close to 
|^|-5/4 ^Yi^Q left of the theoretical shocking distance is the linear evolution regime, where 
AMF increases as the wave accumulates angular momentum. Note that in the linear region, 
a pattern of AMF increase with x should not depend on Mp. The apparent dependence is 
because t in Figure [6] has been scaled by M~^. At /sh the wave AMF stops increasing and 
starts to dissipate. To the right of Ish, in the nonlinear regime the patten of AMF decay 
with X does depend on Mp. However, we remove this dependence and make the AMF decay 
pattern independent of Mp by using t instead of x in Figure El 



5. Effect of Various Numerical Algorithms 

In this section, we investigate the impact of various numerical algorithms on the non- 
linear wave evolution and the subsequent shock formation by looking at their effect on the 
generation of potential vorticity. The algorithms we explore are the same as in paper I, as 
shown in Table (H and values corresponding to our standard case are indicated in boldface. 
Figure [7] shows the radial A( profile for different values of various numerical algorithms, and 
we discuss the effect of each of them in detail below. 



5.1. Solver and its accuracy 

In our simulations we compare two different Riemann solvers — Roe's linearized solver 
jRoelll98lh an d HLLC JTorol Il999l ) — with three different algorithms for the spatial recon- 



struction step ( Stone et al. 20081 ) : second order with limiting in the characteristic variables 



(denoted 2 in this work), which is the predominant choice in literature in this field, and third 
order with limiting in either the characteristic variables (3c), or in the primitive variables 

Table 1. Algorithms of the simulations 



Parameters 



Range 



Riemann solver (flux function) used 

Order of accuracy 

Boundary conditions in y 

Resolution of the simulation (cells per h) 

Planetary potential (see JSj 

Softening length 

Equation of states of the fluid 

Mass of the planet (O.OlMth) 



Roe, HLLC 
2, 3c, 3p 
Outflow, Inflow/Outflow 
64, 128, 256 



(2) 



i(4) ^(6) 



1/8, 1/16, 1/32 
7 = 1, 7 = 5/3 
;.7, 24.2, 11.8, 4.28, 2.09, 1.20, 0.757, 0.369 
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(3p). 

Similar to their effect in the hnear regime (paper I), the two kinds of solvers yield almost 
identical results on the profile. On the other hand, different orders of accuracy do make 
significant differences. 3p accuracy (not shown in this figure) introduces large high frequency 
fluctuations on A( profile compared to 2 and 3c accuracy. For the two spatial reconstruction 
steps with limiting in the characteristic variables, as shown in Figure [7^ 2 accuracy generates 
a much smoother rise of the profile around Igh, compared with the sharp jump in the 
3c accuracy case. In addition, 2 case advances shock formation, causing the numerical Igh 
to disagree with theoretical prediction (|8]). Analogous to Paper I, we find that the effect 
of reducing the accuracy from third order to second order is very similar to reducing the 
resolution by a factor of 2. 

A critical ingredient to any total variation diminishing (TVD) reconstruction scheme, 
as are used in Athena, are the slope limiters used to ensure monotonicity. We have not 
explored the use of different limiter in this work, but this may affect the accuracy of the 3p 
and 3c methods. 



5.2. Resolution 

Resolution has a strong influence on the shock formation, as shown in Figure [7)d. Lower 
resolution considerably accelerates shock formation, causing 1^^ to deviate from the theo- 
retical prediction. Furthermore, the shape of the rising edge of the A^ profile depends on 
resolution. Increasing resolution steepens the edge of the A( curve, and eventually makes it 
almost vertical in the case of 256 /h (which is a sign that the convergence on resolution has 
been achieved). We also find that at low resolution {64/ h in our case), the A^ profile always 
demonstrates a double peaked structure, which is a purely numerical artifact. 



5.3. Planetary potential 

Following paper I, we test the sensitivity of our results to the specific method of potential 
softening, and try three different forms of the softened planetary potential. One of them is 
the second order potential defined as 

= -GM,^^^^„ (23) 

which converges to = —GMp/p at p ^ as (^s/p)^ (which means the fractional error is 
0((rs/p)^) as Tg/p — 7- 0). This is the potential most commonly used in numerical hydrody- 
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namic studies. We also studied the fourth order potential (Eq. [19]), converging to the point 
mass potential as [ts/ p)"^ for p ^ Vg, and the sixth order potential: 



(24) 



converging to = GM^j p as {r^/ p)^ at p ^ rg. We also varied softening length at fixed 
form of the potential. The effect of various $p is shown in panel (c), and that of different 
in panel (d). We find that switching to $p one level higher in accuracy increases A^peak by 
~ 10 — 15%, steepens the edge of the profile, and slightly advances the shock formation. 
Reducing rg by a factor of 2 has quantitatively similar effects. But in general the variation 
of $p and Ts has far less prominent effect on the A^ evolution than the accuracy of solver 
and resolution. 

Based on this discussion we can state that in order to accurately follow the nonlinear 
wave evolution and capture the shock formation, it is crucial to use high order of accuracy 
of the numerical solver, and high spatial resolution. Simulations which do not satisfy these 
criteria may be affected by numerical artifacts and may not be able to resolve properly 
either the nonlinear wave evolution or its consequences — the migration feedback and the 
gap opening. 



6. Nonlinear Evolution in Presence of Explicit Viscosity 

Nonlinear evolution of density waves launched by the low-mass planets can be signif- 
icantly affected by linear damping even if the latter is due to numerical viscosity. This is 
because the time it takes for nonlinearity to have an effect is longer for smaller Mp, which 
enhances the relative contribution of the linear wave damping. 

To explore the effect of linear damping on nonlinear wave evolution, we carry out an 
experiment with explicit Navier-Stokes viscosity. We choose kinematic viscosity v to cor- 
respond to Reynolds number Re = hcs/v = 10^. This is equivalent to having effective 
Shakura-Sunyaev a = v/[hcs) = 10~^. We compare the results with another simulation 
under otherwise identical conditions but without exphcit viscosity (our numerical viscosity 
is at least 10 times smaller in this case, which corresponds to Shakura-Sunyaev a < 10^^). 
The results of the comparison are shown in Figure [H 

As discussed in Paper I, before the shock the numerical AMF and accumulated torque 
calculations should agree with each other, because of the angular momentum conservation. 
After the shock formation, dissipation damps the wave and transfers the angular momentum 
from the wave to the local disk material causing the AMF to drop below the accumulated 
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torque curve. Accumulated numerical torque in our simulations is calculated as: 



X 



Th(x) = / — — dx = dx dySTi—^, (25) 
dx J J oy 





and its asymptotic value in the linear case is (GT80): 

T™ ~ 0.93 (GMp)' (26) 

As shown in the panel (a), while in the inviscid case this deviation happens right at the 
theoretically predicted shocking length (Zsh ~ 4/i), in the viscous case AMF starts to deviate 
from the Th{x) much earlier (at ~ 2.5h). This indicates that dissipation happens earlier in 
the viscous case due to the linear damping, which then must be reflected in the disk density 
redistribution, eventually affecting migration feedback and gap opening. In addition, panel 
(b) shows that in viscous case A( starts to rise from zero much earlier than in the inviscid 
case, also indicating premature dissipation. Furthermore, the rise of in the viscous case 
is very gradual, in contrast to the sharp jump in the inviscid case, and the peak A(^peak is 
also signiflcantly reduced (by a factor of ~ 4). The combination of the two effects makes the 
shock detection very ambiguous in the viscous case. 

We note that normally disks with a = 10"^ are considered to have very low viscosity, if 
not inviscid. However, as we show here, the physics of the nonlinear density wave evolution 
is very subtle so that even an a = 10~^ viscosity can dramatically affect wave damping. 
In simulations which aim at investigating the nonlinear wave evolution, the low spatial 
resolution (< 32 /h, typically employed in the literature) will likely introduce high numerical 
viscosity. This is likely to give rise to inaccurate numerical results in a way similar to our 
viscous simulation shown here, and via the improperly captured back reaction on the disk 
density distribution, affect the description of the migration feedback and gap opening. 



7. Implications for realistic protoplanetary disks 



The analytical and numerical calculations presented in this work directly apply only to 
the case of a 2D, laminar disk, containing a single low mass planet. We now discuss the 
applicability of our results to more realistic protoplanetary disks. 

First, we note that for typical h/vp values of protostellar disks (on the order of 0.1, 
Hartmann et al.lllQQSt IChiang &: Goldreichlll997l ). the size of the radial domain in this work 
(up to lOh on one side) would imply global disk dimensions, which might be seen as con- 
tradictive to the essence of the shearing box simulation. However, our numerical calculation 
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is an idealization needed for che cking th e GROl theory. Extensions of the GROl theory to 
the global case are available in iRafikovl ( l2002a( ). with which global simulations should be 
compared. Here we just point out that qualitatively things are the same in the global case. 

Excitation of density waves by planets in fully three-dimensional ( 3D) disks has been pre- 
viously investigated by a riumbe r of authors ( iLubow fc Pringlelll993l : iKorycansky fc Pringle 



1995 



Takeuchi fc Miyamal 119991 ). These studies have generally shown that disks with verti- 



cal therr nal stratification do not support the modes sirnilar to the modes existing in a purely 
2D disk jKorvcanskv fc PrindelEggsl : Ihubow fc 0gilvielE998l : l0gilvie fc Lubowlll999[ ). Thus, 
the results of our study cannot be directly applied to such thermally stratified disks. 

However, protoplanetary disks at a dvanced stages of planet formation are expected 
to be passive (jChiang &: Goldreichl Il997l ). heated predominantly by the radiation of their 
central stars with only negligible contribution from accretional energy release. Because of the 
extern al illumination such disks are expected to be verti cally isothermal. It was previously 
found (ILubow fc Pringldll993l : iTakeuchi fc Miyamalll999l ) that in vertically isothermal disks 
planetary gravity always very effectively excites the two-dimensional mode with no vertical 
motion, which is similar to the density wave in a 2D disk studied here. Even though other 
modes with non-zero vertical velocity perturbation are also excited by the planet in 3D 
isothermal disks, these modes are found to carry only a small fraction of the total angular 
momentum flux (ITakeuchi fc Miyamal Il999l ). Thus, density wave excitation and nonlinear 
dissipation in realistic passive protoplanetary disks should be very similar to the picture 
outlined in GROl and this work. 

The assumption of the laminar inviscid background flow in the disk used in this work 
may be violated in realistic protoplanetary disk if it is strongly turbulent, e.g. due to the op- 
eration of the magnetorotational instability (MRI). However, this instability and associated 
deviations from the purely laminar flow are expected t o be greatly suppressed in the so-called 
"dead zones" of protoplanetary disks (jGammielll996l ). which are so weakl y ionized that the 
opera tion of MRI is not supported there because of non-ideal MHD effects (IFleming &: Stone 
20031 ) . These zones are expect ed to extend over th e several- AU wide region where planet 
formation is expected to occur (ITurner fc Sandl2010l ) and in this part of the disk our results 
obtained under the assumption of laminar flow should be applicable. One should also men- 
tion that because of the suppressed viscous dissipation inside the dead zone the assumption 
of a vertically isothermal thermal structure of these regions is well justified. 

Protoplanetary disks may easily host a number of protoplanets simultaneously and one 
may wonder whether the overlap of the density waves excited by different objects can have 
some effect on their evolution. This effect is not important in the limit of low mass planets 
Mp < Mth explored in this work since the density waves excited by such objects are only 
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weakly nonlinear. As a result, the possible overlap of such density waves is, to a good approx- 
imation, just a linear combination of the independent density waves excited by individual 
planets. To summarize, the results of our work should directly apply at the very least to 
the density wave evolution in low viscosity regions ("dead zones") of passive protoplanetary 
disks heated predominantly by their central stars. 

8. Sumniciry 

We present a numerical study of the nonlinear density wave evolution, based on 2D 
local shearing box hydrodynamical simulations using the grid-based code Athena. The non- 
linear wave evolution plays an important role in the disk-planet interaction and planetary 
migration. It provides an efficient and robust wave damping mechanism, working in regimes 
where linear dissipation mechanisms fail, and leading to gap opening by very low mass plan- 
ets (~ M®). Nonlinear wave evolution leads to shock formation, wave damping, and transfer 
of the angular momentum carried by the wave to the local disk material. The latter drives 
global disk evolution by redistributing the mass in the disk, and leads to the migration 
feedback and gap opening. 

We analytically study the potential vorticity generation at the shock location, and verify 
our calculation numerically. The scaling relation between A^^peak and Mp derived from the 
simulations is very close to the analytical result. We use the jump in potential vorticity as 
a way of pinpointing the shock location, based on which we numerically confirm with high 
accuracy (better than 10% in most of the Mp range) the theoretically predicted /gh — Mp 
relation for Mp varying by 2.5 orders of magnitude, from few lunar masses to several Earth 
masses at 1 AU. In addition, the theoretical dependence of the shocking length l^h — Mp 
relation on the equation of state of the gas is also verified. 

We investigate the evolution of the density wave profile in the post-shock regime, and 
observe its evolution towards the N-wave shape. We verify the theoretical prediction for the 
evolution of the peak amplitude and the azimuthal width of the N-wavc. The post-shock 
decay of the angular momentum fiux carried by the wave also agrees with theory. 

Furthermore, the effect of various numerical algorithms on the simulation results is 
explored, including numerical solver and its accuracy, resolution, planetary potential and 
the softening length. We find that in order to accurately follow the nonlinear wave evolution 
and capture the shock formation, high order of accuracy for the solver and high resolution 
are crucial. Simulations which do not satisfy these criteria will have trouble resolving shock 
formation and post-shock wave evolution driven by nonlinearity. The latter may further 
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affect ttie consequences of tfie nonlinear wave evolution, such as the migration feedback and 
gap opening by low mass planets. 

In addition, we find that the linear viscous damping strongly affects shock formation and 
the nonlinear wave evolution. Using an experiment with explicit viscosity, we find that the 
viscosity at the level of a ~ 10~^ can strongly modify the generation of potential vorticity at 
the shock, and accelerate the dissipation of the angular momentum carried by the wave. In 
low resolution simulations, the resulting high numerical viscosity may lead to similar effects 
on the nonlinear wave evolution and negatively impact the results on disk-planet interaction 
and planetary migration. 
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Fig. 1. — Left: Density profiles across tlie density wake for the Mp = 2.09 x 10~^Mth 
(corresponding to Igh ~ 4/i) case at two pre-shock locations (1.33/i and 2.67h, thin solid 
curves), four post-shock locations (5.33/i, 6.67/i, 8.0h, and 9.33/i, dashed curves), and the 
theoretically predicted shocking length {4h, the thick solid curve). Right: Scaling of A and 
w with coordinate t{x) defined by Eq. ([5]) in the post shock region. Theoretical scaling 
relations (Eq. [TUi) are over-plotted with arbitrary normalization. 
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Fig. 2. — A typical snapshot of the potential vorticity perturbation A<^ (normalized by 
n/Eo) in shearing sheet coordinates in our simulation (only half of the simulation domain 
is shown). A planet with Mp = 4.28 x 10~^Mth (corresponding to ~ <ih) is located at 
X = y = 0. Fluid enters from the upper right and the lower left boundaries. Note the 
vorticity generation at the shock position, with the amplitude of A( decaying far from the 
planet (in x). Vorticity perturbation in the horseshoe region is not related to the shock. 
Narrow transient features in the preshock region are due to fluid activity in the horseshoe 
region. 
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Fig. 3. — Panel (a): Numerical peak A^peak amplitude as a function of Mp. Best fit relation 
C oc Mp-^^ and the theoretical scaling relation A(peak oc M^-^ (Eq. [17]) are shown. Panel (b): 
Radial profile of A( (scaled by i^/So) as a function of x (scaled by Ish)- Different curves for 
different Mp have been scaled by (Mp/Mth)^'^^ (the numerically best fit scaling factor) to 
remove the dependence on Mp. See the discussion in §4.21 




Fig. 4. — Radial profile for two simulations with otherwise identical parameters (including 
the adiabatic sound speed) but different EOS (7 = 1 and 7 = 5/3). Mp = 2.09 x IQ-^Mth 
in both cases, which corresponds to Igh ~ 4/i for 7 = 1 and Ish ~ 3.6h for 7 = 5/3. 
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Fig. 5. — Numerical Igh as a function of Mp (dots) as well as the theoretical prediction of 
GROl (Eq. ([8]), solid line). Ish is determined as the midpoint of the jump at the shock, 
as we discussed in § 14.21 For the five low mass cases we use resolution=256//i and = h/32, 
and for the three high mass cases we use resolution=128//i and rg = h/16. The smallest and 
largest Mp here correspond to a few Lunar and Earth mass at 1 AU. 
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Fig. 6. — Numerical result for the AMF decay after the shock formation for three Mp. The 
two axes are scaled to facihtate direct comparison with Fig. 3 in GROl (see §4.51 for details). 
For the two low mass cases we use resolution 256/ h and = h/32, and for the highest mass 
case we use resolution 128 /h and rg = h/16. Theoretical asymptotic AMF decay scaling 
relation (Eq. (fTTj) ) at x l^h from GROl and the position of x = l^h are indicated by dotted 
lines. 
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Fig. 7. — Radial profile of for simulations with different order of accuracy (a, with 
resolution=128//i, = h/16), resolution (b, rg = h/16), $p (c, with resolution=128//i, and 
Tg = h/16), and Vg (d). Other numerical algorithms which are not mentioned are drawn from 
our fiducial choices (©. For all simulations we use Mp = 2.09 x 10~^Mth (corresponding to 
^sh ~4/i). 
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Fig. 8. — The effect of explicit viscosity in our runs. The viscous simulation is done with 
an explicit Navier-Stokes viscosity (Shakura-Sunyaev a = u/lhcg) = 10~^), and the inviscid 
run is our standard simulation without explicit viscosity. Both simulations arc done with 
Mp = 2.09 X 10~^Mth (corresponding to /gh ~ 4/i). Panel (a) shows the numerical AMF and 
torque calculation, and panel (b) shows the radial A( profile. 



